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For a bipartite honeycomb lattice, we show that the Berry phase depends not only on the shape ol 
the system but also on the hopping couplings. Using the entanglement entropy spectra obtained by 
diagonalizing the block Green's function matrices, the maximal entangled state with the eigenvalue 
A m = 1/2 of the reduced density matrix is shown to have one-to-one correspondence to the zero 
energy states of the lattice with open boundaries, which depends on the Berry phase. For the 
systems with finite bearded edges along x-direction we find critical hopping couplings: the maximal 
entangled states (zero-energy states) appear pair by pair if one increases the hopping coupling h 
over the critical couplings h c s. 
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Experimental studies on graphenefH and cold atomsQ 
have inspired new interests in quantum properties of elec- 
trons and atoms on the honeycomb lattice. Due to its pe- 
culiar energy dispersion, which contains two Dirac points 
as Fermi surfaces, the honeycomb lattice offers particu- 
lar physical phenomena which can not be observed on the 
square lattice. In particular the edge states with zero en- 
ergy may appear when special boundaries are present. 
The existence of the edges state relies on the shape of 
boundaries, that is, for zigzag and bearded edges there 
exist zero-energy edge states, on the contrary the arm- 
chair edges offer no zero-energy edge modes according to 
their Berry phase of the occupied band Q . Therefore the 
topological structure of the system strongly influences the 
existence of zero-energy edge modes. 

On the other hand, recent development on quantum 
information applied in condensed matter theory provides 
a new tool to investigate quantum phenomena [4] . Espe- 
cially the entanglement entropy: Sa = — Tr pa log 2 pa > 
where pa — Tr# \^>ab}(^ab\, has been widely used to 
measure the bipartite entanglement for a given pure state 
(ground state) I^ab) of a bipartite AB consisting of sys- 
tem A and environment B. In the last several years this 
quantity has been successfully used to interpret the na- 
ture of the quantum criticality [B| due to the fact that the 
entanglement entropy diverges at the quantum critical 
points of the second order phase transition. The scaling 
law of the entanglement entropy even provides a rule to 
tell different quantum phases [6|. Furthermore, the en- 
tanglement entropy spectrum for free fermions, defined 

as the set of A m , where pa = ® m 'V 1 , can 

U 1 A m 

be used as a tool to investigate physical properties such 
as disorder lines [7[, Berry phase [8j and zero-energy edge 
states [8j. While the entanglement entropy Sa = Yl m 
where 

S m = -A m log 2 A m - (1 - A TO ) log 2 (1 - A m ), (1) 





FIG. 1: (Color online) a honeycomb lattice with different 
shapes for system A. (/U r -green rectangle) two bearded 
edges, (A zz - red rectangle) two zigzag edges, (A ar -blue rect- 
angle) two armchair edges. The hopping constant along x- 
direction is h, along y-direction 2. On the right hand side 
the transformation from a honeycomb lattice to a brick-type 
lattice is shown. 



contains more global information about the system, the 
entanglement entropy spectra offer a new method to ob- 
serve the microscopic quantum phenomena. 

In this Letter we explore the bipartite entanglement 
entropy spectra of a tight-binding honeycomb lattice for 
different shapes of the system A, as shown in Fig. [T] We 
show that not only the shape of the edges but also the 
hopping coupling h along the x-direction influences the 
Berry phase x, which value decides if the system has max- 
imal entangled states with A m = 1/2 or not [8]. In the 
previous literature, only the edges with infinite size have 
been studied Q , while in this Letter we analytically and 
numerically study the edges with finite sizes. For bearded 
edges, the maximal entagled states appear two-by-two 
only if the hopping coupling h passes by critical hopping 
couplings h c s. These results are verified numerically us- 
ing the entanglement entropy spectra by diagonalizing 
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FIG. 2: (Color online) Entanglement entropy spectra (red 
circles: R z = 0, blue dots:.R z = 0.01) (a) h = 2, L x = 3 
and L y — 200 (two bearded edges); (b)h = 2, L x = 201 and 
L y = 4 (two armchair edges); (c) h = 2, L x = 11 and L y — 
200 (two bearded edges); (d) h = i,L x = 11 and L H = 200 
(two bearded edges). The insets show the number of maximal 
entangled states. 



the block Green's function matrix. 

We consider a honeycomb lattice with anisotropic hop- 
ping constants as shown in Fig. [TJ The hopping coupling 
along the y direction is defined as 2 of magnitude, while 
the other two couplings are replaced with a general value 
h. Without changing the lattice topology, a honeycomb 
lattice can be transformed into a brick-type lattice^, 
therefore the Hamiltonian can be rewritten as 



E 

x,y 



{hcl. y c x+l , y +[l + {-l) x + y ] cl y c x<y+1 + h.c.}. 

(2) 

The dispersion relation for such Hamiltonian is: e(k) = 
2(1 + h 2 cos 2 k x + 2hcosk x cosfc^) 1 / 2 . For < h < 1, the 
system is gapped, while for ft > 1 we find two Dirac 
Fermi points located at (k x ,k y ) = (± cos" 1 (— 1/h), 0). 
For h = 1, the two Dirac zero-energy modes merge into 
a confluent point (k x , k y ) = (ir,0). The zero-energy edge 
states for bearded edges appear only if the bulk has two 
Dirac points. The reason is described as follows. 

The edges on the honeycomb lattice can be produced 
by imposing different open boundaries: bearded, zigzag 
or armchair edges. To further see the influence of the 
boundaries, the sub-lattices 7 and 5 can be first labeled 
as o and • in Fig. [TJ Defining a pair of annihilation 
operators = (c*y,c$) x and Fourier transforming the 
Hamiltonian ©, we can rewrite the Hamiltonian @ in 
the momentum space as 



n = -J2 4[R(k)-<7] Ck , 



(3) 



keflz 



where a = (cr x ,a y ,a z ) as Pauli matrices and R(k) = 
(R x , Ry,Rz) € R 3 . If we consider systems without onsite 



potentials, one can always rotate R(k) to lie on a two- 
dimensional plane by applying a global 5*0(3) unitary 
transformation, thus one can define R(k) = (R x ,R y ). 
In this parametrization, the energy e(k) = ±|R|. The 
form of R(k) is given by the choice of the unit cell 
in Fourier transformation along the edges Q. (a) For 
a bearded edge along x-direction, R(k) = (h[cosk y + 
cos(k x — k y )}+2, h[smk y — sin (k x — k y )]); (b) for a zigzag 
edge along x-direction, R(k) = (2cos(fc;E — k y ) + h[l + 
cos k x ] , h sin k x — 2 sin (k x — ky)); (c) for a armchair edge 
along y direction, R(k) = (h[cosk x + cos(k x + ky)] + 
2, h [sin k x + sin (k. x + k y )]). The choice of the unit cell in 
Fourier transformation is not unique, however, different 
choices lead to the same topology for the edge states ■ 
By fixing a wave factor parallel to the edges and investi- 
gating the loop of R as a parameter of the perpendicular 
wave vector changing from — ir to n, the topology of the 
system can be obtained. In case that the loop £ of R con- 
tains the origin O in the R space, the Berry phase (or 
Zak's phase) x, defined as a line integral of the curvature 
of the filled band, is ir. In this case, due to the fact that 
one can continuously deform I into a unit circle without 
crossing the origin, the topological argument ensures us 
that the original Hamiltonian corresponding to I contains 
at least one zero-energy edge state 3|. Therefore for (a) 
bearded edges along x-direction, where k x is fixed and k y 
serves as parameter of the loops, there exists a zero en- 
ergy state only if — 2cos _1 (1/h) < k x < 2COS" 1 (1/h). 
Similar situation happens for (b) zig-zag edges. The 
zero-energy edges states appear under the condition that 
-7T < k x < -cos" 1 (2/h 2 - 1) or cos" 1 (2/h 2 - 1) < 
k x < 7T. On the other hand, for (c) armchair edges, 
no zero-energy states exist due to the fact that no loops 
will encircle O with fixed k y by changing k x . Therefore 
for the zigzag and bearded edges the hoping coupling h 
changes the range of k x where the Berry phase \ — ir. 

The consideration for entanglement entropy is some- 
how different: we no longer cut the system with an 
open boundary, but partition that into two parts: sys- 
tem A and environment B. By using von Neumann en- 
tropy Sa, one can figure out how the system A entan- 
gles with the environment B. However, there exists a 
one-to-one correspondence between the zero-energy state 
for the Hamiltonian with edges and the maximal entan- 
gled state for the bipartite system. The reason is as 
follows. Consider the whole bipartite AB consisting of 
N sites (or modes), with n sites of system A and no- 
tice that in our calculation only thermodynamic limit 
(iV — > 00) will be taken. The reduced density matrix 
Pa = Tr^po where po = |^'ab)( v I / j 4b| can be obtained 
by determining the matrix elements of the full density 
matrix with respect to coherent states and integrating 
out the variables of the environment B [7[ and can be re- 
lated to the eigenvalues A m of the block Green's function 
matrix[ll[ Gys^i - r,) = Trpc( ri)7 )c| r . )(5) , where (1-1,7) 
and (tj,6) belong to system A. The entanglement en- 
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FIG. 3: Position labeling for wavefunctions of the honeycomb 
lattice. (a)L x — 3 (b) L x — 11. 



tropy therefore takes the form as Eq. ([T]). In this case, 
the Green's function matrix can be Fourier transformed 
as G 75 (r, - r,-) = A^ 1 £ keBZ e-^'-^G^k), where 
G 7 5(k) is calculated as 



GV(k) = 2 



R -a 

R 



(4) 
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We regard G^si^i — as an effective Hamiltonian for 
obtaining the entanglement entropy. In the case of tak- 
ing the whole infinite plane with periodic boundary, i.e. 
no sites in the environment B, the G Hamiltonian has 
the same set of eigenfunctions as the original Hamilto- 
nian, however, the eigenvalues are either 1 or 0, which 
means 5 = 0, the system is not entangled at all. On the 
other hand, while B is not empty, the nontrivial bound- 
ary states appear, which eigenvalues are not 1 or 0. This 
boundary states, which eigenvalues < A m < 1, cause 
the nonzero entanglement entropy, therefore we can say 
that those are the most important states for calculating 
the entanglement. 

Comparing the Hamiltonian ^ with G 7 a(k) they 
almost take the same form except a constant and a pos- 
itive normalization factor R, therefore we can conclude 
that they should share the same topology, that is, if R 
encloses origin in the parameter space which causes Berry 
phase equal to ir, the zero-energy state appears for the 
Hamiltonian , while for the block Green's function ma- 
trix we obtain one special eigenvalue A m = 1/2 [8]. They 
also share the same eigenfunctions with different eigen- 
values. We call the state with A m = 1/2 a maximal en- 
tangled state due to the fact that the resulting S m = 1. 
According to the discussion above, only the reduced den- 
sity matrix of the system A with bearded and zigzag 
edges has maximal entangled eigenstates, while for the 
armchair edges there exist no maximal entangled states. 

In order to investigate this effect we numerically diago- 
nalize the block Green's function matrix G 7 s(^i — ry) for 
the system A with finite size, as shown in Fig. [TJ which 



L x and L y are defined as follows. Transforming the hon- 
eycomb lattice into the brick-type lattice and projecting 
the system sites into x and y direction, L x (L y ) is de- 
fined as the total number of the projected sites along 
the x(y) direction. The real noninteracting edge state 
only appears in the thermodynamic limit, i.e. L x L y 
or vice versa. For example the partition of system A\, r 
and A zz m Fig.jTJhave either two bearded (At r ) or zigzag 
edges (A zz ). On the other hand there exist two armchair 
edges for the system A ar . Here we only consider the edges 
far apart from each other so that the two edges do not 
interact with each other (independent edges). Therefore 
we expect that for the partitions Af, r or A zz there ex- 
ist maximal entangled states (A m = 1/2), while for A ar 
the maximal entangled states do not appear. (Due to 
the fact that bearded and zigzag edges have the same 
properties, we will only calculate the entanglement spec- 
tra for the bearded edges in the following.) The upper 
panels of Fig. [5] shows the entanglement spectra A m of 
a system size: (a) L x = 3, L y = 200 (Ab r -type) and (b) 
L x = 201, L y — 4 (j4 Qr -type) for a graphene (h = 2). 
For the system with two independent bearded edges, we 
obtain two maximal entangled states with A m = 1/2 (see 
the inset), while there are no maximal entangles states 
for the system L. x 3> L y (on the right hand side of Fig. [JJ 
due to the fact that ^4 ar -type system has two indepen- 
dent armchair edges. Fig. [5] (c) shows the spectra for 
L x = 11, L y = 200 and h = 2. From the inset we ob- 
serve 6 eigenvalues close to 1/2. In Fig. HJd) the system 
size is the same as (c) but the hopping coupling h = 4. 
We obtain 10 maximal entangled states (Due to finite- 
size effect, two eigenvalues of the states are not precise 
1/2, however we still count them as maximal entangled 
states). The number of maximal entangled states not 
only depend on the shape of the system but also on the 
hopping coupling h. We will discuss that as follows. 

For general h, the wave function of the maximal en- 
tangled state is the same as the zero-energy mode. The 
zero-energy wavefunction is obtained by the tight-binding 
Schrodinger equation (SE) — y\ tjipj = Etpi — 0, where 
j are neighbor sites of i, tj — h for the vertical hop- 
ping couplings and tj — 2 for the longitudinal hopping 
constants in the brick-type representation. Fig. [3] (a) 
shows the structure of the system with L y 3> L x = 3. 
To obey SE, a solution with all values of the wavefunc- 
tion for o sites in Fig. [3] equal to zero is expected. For 
each partition symmetric and antisymmetric wavefunc- 
tions with respect to the middle line of the system in 
the x direction can be separated to two different groups. 
For symmetric solutions, we set that [■0Ai,V , a 2 ] = [<A) ( / ) ]- 
According to SE, i/jb 1 = —(2/h)tp and [^pcni'c^] — 
[(2/h 2 )t/j, (2/h 2 )<j)]. Therefore there only exist maximal 
entangled states (zero-energy states) if 2/h 2 < 1, other- 
wise the wavefunction will explode with the increasing 
L y . We obtain a critical h c = y/2: the maximal entan- 
gled states appear only if h < h c . On the other hand 
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FIG. 4: (Color online) (The main figure) Numerical result of 
the nearest eigenvalues close to 1/2 with varying h for L x = 11 
and L y — 200. (Lower inset) Zooming for A m near 1/2. There 
exist five critical couplings h c . (Upper inset) Analytical result 
for h c with different sizes: L x = 3, 5, 7, 9, 11. 

the antisymmetric solutions [tpAu^pAal = f an f° r 

L x = 3 due to the fact that the value on site B\ directly 
sets (j> = 0. 

Another example we study here is L x — 11, as shown 
in Fig. [3] (b). For the symmetric wavefunctions, we 
can first set the values [V>aJ = [1, otti a 2, «2, cxi, ^}4> and 
[ipBi] — 02, fa, fa, fa]<t>- Due to the fact that the 
structure of the layer C\ is the same as Ai, the val- 
ues of the wavefunction can be set as [V>cJ = x [^Ai]- 
To obtain the critical couplings h c s we have to know 
all these values an , fa and x by using SE in different 
o sites. The SE for the site between A\ and B\ gives 
fa = —2/h. The site among A 2 ,B\ and B-i leads to 
the equation fa = 2(1 — ax)/ h. In the same manner 
we obtain fa = 2{a 2 — ot\ — l)/h. Using these relations 
and the SE for the o site between B2 and C 2 , the first 
value for x is found: x = 4/[(l + ax)h ]. From this re- 
lation the critical coupling is obtained by setting x = 1, 
i.e. h c = 2/ + a\. h c is therefore obtained if u\ is 
known. On the other hand, different relations for x can 
be obtained by subsequently calculating the SE for the 
o sites among Bi and 6^: x = 4(ai — l)/[(ai + a 2 )h 2 ] 
and x = 2(1 + a 2 — ai)/(a 2 h 2 ). These three x should 
be identical, we first obtain a 2 = a\ — a\ — 1 and finally 
the polynomials for a\ is found : a\ — 3a 2 + 2 = 0. This 
indicates that a\ — 1 or ai — 1± \/3. For the sym- 
metric eigenfunctions we obtain three critical couplings 
h c = \/2, 2 V2- \/3 and 2 v2 + \/3 corresponding to dif- 
ferent ax- 

The antisymmetric eigenfunctions can be also found 
in a similar way. We can first assume that [ipAi] = 
[l,ai,a 2 ,-0!2,-a!i,- !]<?!>, [ip Bi ] = [fa, fa, 0, -fa, -fa]<j> 



and [ipd] = x[ipAi]- From the layer between A, and Bi 
one obtains fa = —2/h and fa = 2(1 — a.2)Jh. The 
SE for the site among A3, B2 and B3 gives the rela- 
tion as a2 = (a± — 1). Two relations for x can be ob- 
tained by solving the SE for the sites between Bi and GV 
x = 4/[(l + ax)h 2 ] and x = 4(ai - l)/[(ai + a 2 )h 2 ]. 
Using these relations we finally obtain a\ — 0, 2 and 
h c = 2y/3/3 and 2. Fig. H shows the some eigenvalues 
A m nearest 0.5 of the system L x — 11 and L y — 200 with 
bearded edges along x-direction as a function of h by di- 
agonalizing the block Green's function matrix. The lower 
inset shows that the system has 5 different h c s and they 
coincide with the analytical values we discussed above. 

The upper inset in Fig. 2] shows the analytical values 
of h c s with varied L x . Firstly, for a fix h, the eigenvalues 
increase as a function of L x . In the case of L x — > 00, 
those eigenvalues form a continuum, i.e. form a flat band 
of zero surface states, as shown in [9( for graphenes (h = 
2). The number of eigenvalues can be approximated as 
N sa 2/ttL x cos -1 (1/h) for large L x . This accounts for 
the increasing number of the eigenvalues as h increases. 

Secondly we found that there are some repeated h c s 
with growing L x . For example: for L x — 3,7,11, the 
system has the same h c = y2- The reason is that one can 
build either a symmetric or antisymmetric wavefunction 
for L x — 7 and 11 using a fundamental block L x — 3 
similar to the toy Lego : a positive and a negative block 
of L x — 3 can build a system with L x — 7; a sandwich- 
like system with two L x — 3 positive blocks on both 
sides and a negative block in the middle has exact the 
size L x = 11. Similarly, for L x — 5 and L x = 11 there 
exists the same h c = 2. It is clear that a positive and 
a negative block with L x — 5 can build a system with 
L x = 11. 

An interesting question arises: if considering a sys- 
tem with onsite potentials, i.e. if putting R = 
(R x ,R y ,R z ),R z = A, do the edges states still robustly 
exist? By adding the onsite potentials, the chiral symme- 
try is destroyed, there exist no maximal entangled states. 
The blue lines and dots of the inset in FigJ5] shows the 
small deviation of A = 1/2 by putting A = 0.01. How- 
ever, due to the fact that the deviation is small, we can 
conclude that the states with eigenvalues closed to 1/2 
are still edge states. That means, the edge states robustly 
exist even with broken chiral symmetry. 

In the summary, we have used the entanglement spec- 
tra to study the edge states of a honeycomb lattice. The 
dependence between the Berry phase and hopping cou- 
plings has been found. The one-to-one correspondence 
between the maximal entangled state for a bipartite sys- 
tem and the zero energy state has been proved numer- 
ically by diagonalizing the block Green's function ma- 
trix. We found there exist critical coupling couplings h c s 
for the system with finite edges: the maximal entangled 
states pair-by-pair increases every time when the h jumps 
over h c s. 
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